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ABSTRACT 

A key goal of many Cosmic Microwave Background experiments is the detection of 
gravitational waves, through their B-mode polarization signal at large scales. To ex- 
tract such a signal requires modelling contamination from the Galaxy. Using the Planck 
experiment as an example, we investigate the impact of incorrectly modelling fore- 
grounds on estimates of the polarized CMB, quantified by the bias in tensor-to-scalar 
ratio r, and optical depth r. We use a Bayesian parameter estimation method to 
estimate the CMB, synchrotron, and thermal dust components from simulated ob- 
servations spanning 30-353 GHz, starting from a model that fits the simulated data, 
returning r < 0.03 at 95% confidence for an r = model, and r = 0.09 ± 0.03 for an 
r = 0.1 model. We then introduce a set of mismatches between the simulated data 
and assumed model. Including a curvature of the synchrotron spectral index with 
frequency, but assuming a power-law model, can bias r high by '^ la (Sr ^ 0.03). 
A similar bias is seen for thermal dust with a modified black-body frequency depen- 
dence, incorrectly modelled as a power-law. If too much freedom is allowed in the 
model, for example fitting for spectral indices in 3 degree pixels over the sky with 
physically reasonable priors, we find r can be biased up to ~ 3a high by effectively 
setting the indices to the wrong values. Increasing the signal-to-noise ratio by reducing 
parameters, or adding additional foreground data, reduces the bias. We also find that 
neglecting a '^ 1% polarized free- free or spinning dust component has a negligible ef- 
fect on r. These tests highlight the importance of modelling the foregrounds in a way 
that allows for sufficient complexity, while minimizing the number of free parameters. 



1 INTRODUCTION 

Extraction of the polarized Cosmic Microwave Background 
signal at large scales is hampered by significant levels of 
polarized Galactic emission. The two dominant components 
are synchrotron and thermal dust, polariz ed due to the co- 
herent magnetic field in the Galaxy fe.g.. |Page et al.ll2007l : 
[Fraisse et al. 2008). For an experiment observing at mul- 
tiple frequencies, one method of separating the signals is 
to parameterize the synchrotron and dust, and to fit for 
these components, in addition to the CMB, over the re- 
gion of sky where Galactic emission is lowest. While demon- 
strated to work for E-mode polarizati on (Page et al. 2007, ; 
iDunklev et all l2009al : [Gold et al.ll2009l '). the signal of inter- 
est is the much smaller B-mode signal from inflation (e.g.. 



armitage-caplan@physics.ox.ac.uk 



iBasko fc Polnarev |[T980l : iBond fc Efstathiou II 1984 ). A con- 
cern with using such methods is that an incorrect model can 
lead to bias in the estimated CMB signal. 

The Planck satellite mission, launched in May 
2009, is measuring the polarization signal of the 
CMB in seven channels over the f r equency range 30- 
353 GHz JPlanck CoUaborationl l2006l : iTauber et aiTl2010l : 
iPIanck Collaboration 11120111 ). While Planck will produce 
polarization data, which offer a multitude of opportunities 
including possible recovery of inflationary B-modes at large 
scales and greater understanding of the polarized nature of 
Galactic foregrounds, it also comes with great challenges. 
For an all-sky experiment like Planck, component separa- 
tion of the polarization signal is more difficult than for the 
temperature counterpart, in part because the ratio of the 
foreground signal to CMB signal is higher. 

In many simulated tests of component separation, the 
simulations of the Galactic emission are well matched to 
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the model used to describe them. Using a Bayesian compo- 
nent separation method which allows us to assume different 
models of the Galactic signal, we explore the effect on the 
recovery of the CMB in varying scenarios of mismatch be- 
tween the model and simulation. We use the recovered CMB 
map and its covariance to estimate two cosmological param- 
eters: the optical depth to reionization, r, and the tensor-to- 
scalar ratio, r. In this way, we can directly quantify the bias 
generated in the parameter estimation as a result of any 
particular model-simulation mismatch. Both the Bayesian 
component estimation method and the simulated skies used 
in t his paper were first used and d escribed in a previous pa- 
per, [Ai7nitage:£a2lan_eLalJ (|201ll ). There we examined the 
prospects for large-scale polarized map and cosmological pa- 
rameter estimation with simulated Planck data for a single 
model-simulation combination. This paper is a natural ex- 
tension in which we use the same methods to recover maps 
and estimate parameters, while varying the simulated data 
and separation model. 

In ^ we provide a brief overview of the Gibbs sampling 
method, the subsequent processing of the sampled distribu- 
tion, and the likelihood estimation method. In iJH we explain 
how the data are simulated. A detailed account of the mis- 
match tests that we examine, and their resulting parameter 
estimates, is then presented in iJJ] We then discuss the re- 
sults and methods for mitigating possible biases in ij5] and 
conclude in ^ 



2 METHOD 

In the Bayesian parameter estimation method of foreground 
removal, the emission models of the CMB and foregrounds 
are parametrized based on our understanding of their fre- 
quency dependence. Focusing on polarization analysis, a 
sampling method is then used to estimate the marginal- 
ized CMB Q and U Stokes vector maps (and additionally 
the marginalized foreground maps) in every pixel over the 
sky. In general this extends template-removal methods to al- 
low for spatial variation of the foreground spectral indices, 
and was first used to clean WMAP polarization data in 
iDunklev et al . (2009a). In this analysis, we use HEALPix 
i Gorski et al.l:2005. ') A'^ido = 16 maps containi ng N^ = 3072 
pixels . We use a code cal l ed Co mmander (see lEriksen et al.l 
l|2006t) and lEriksen et al.l (|2008l ')') to perform the Gibbs sam- 



pling. The sampled distribution is then processed into a 
mean map and covariance matrix. Finally, we perform a 
likelihood estimation for the two cosmological parameters, 
r and r. 



2.1 Bayesian Estimation of sky maps 

By Bayes' theorem, the posterior distribution for parame- 
ters, s, given a set of maps, d, can be written as 



P(s|d) oc P(d|s)P(s) 



(1) 



with a prior distribution for the model parameters, P(s). 
The Gaussian likelihood of the observed maps is given by 

- 21nP(djs) = ^[d. - s,.]^N;'[d, - s.] (2) 



where di, is the observed sky map at frequency v, and N^ 

is its covar iance matrix. 

As in iDunklev et al.l (l2009al '): lArmitage-Caplan et al.l 

|201l| ). we assume that the polarized Galactic emission 
is dominated by synchrotron and dust emission, arising 
due to the orien tation of the Galactic magnetic field (e.g., 
iPage et al.ll2007l ). We define a parametric model for the to- 
tal sky signal in antenna temperature for a three-component 
model (fc = 1 for CMB, k — 2 for synchrotron emission, and 
fc = 3 for thermal dust emission) as 



s^ = ai(i^)Ai + cy.2{iy;l32)A2 + a-i{iy;fy)A3 



(3) 



where Ak are amplitude vectors of length 2Np and Oik{i^; Pk) 
are diagonal coefficient matrices of side 2Np at each fre- 
quency. 

Once our model, and priors on the model parameters, 
are defined, we estimate the joint CMB-foreground posterior 
P(A,/3|d) from which we can then obtain the marginalized 
distribution for the CMB map vector. 



(4) 



p(Ai,d)= / p(A,/31d)dA2dA3d/3 



and similarly for the other model parameters. 

For the multivariate problem that we are considering 
Gibbs sampling draws from the joint distribution by sam- 
pling each parameter conditionally as follows 

^.+1 



P(A|/3,d) 
P(/3|A,d). 



(5) 
(6) 



We use Commander to implement the sampling of the 
amplitude-type and spectral index parameters. Comman- 
der is a flexible code for joint component separation and 
CMB power spectrum estima tion; the reader is directed to 
lArmitage-Caplan et al.l (|201lh for a full description of its use 
for sampling only the sky signal. 



2.2 Likelihood estimation of cosmological 
parameters 

The product of a Bayesian parametric map estimation 
method is both a CMB map (which is taken to be the 
mean map calculated from the Gibbs chain after some burn- 
in) and a covariance matrix (which can be estimated from 
the marginalized posterior distribution) and together these 
products can be used to place constraints on cosmological 
parameters. We compute the likelihood of the estimated 
maps, given a theoretical angular power spect rum, using 
the exact pixel-likelihood method de scribed in iPage et al.l 
(|200'iD : lArmitage-Caplan et all l|201ll ). 

The two cosmological parameters constrained by the 
large scale CMB polarization signal are the optical depth to 
reionization, r, and the tensor-to-scalar ratio, r. The signa- 
ture of reionization is at £ ^ 20 in Cf^ where the amplitude 
of the reionization signal is proportional to r^. The tensor- 
to-scalar ratio r directly scales the Cf^ power spectrum 
and is best probed at two angular scales: at the low £ ^ 20 
'reionization bump' before Cf^ due to lensing dominates, 
or at the smaller scale £ ~ 100 'recombination bump' where 
foregrounds are expected to be lower but lensing is a con- 
taminant. In this study we are considering constraints from 
the large-scale reionization bump, using ~75% of the sky. 



By varying only the optical depth to reionization, and 
fixing the temperature anisotropy power at the first acous- 
tic peak {£ = 220) , we calculate the likelihood for each value 
of T. Separately, we vary only the tensor-to-scalar ratio, and 
calculate the likelihood at each value of r. The resulting one- 
dimensional distributions for r and r then include marginal- 
ization over foreground uncertainty. To account for imper- 
fect foreground cleaning in the Galactic plane, we apply a 
Galactic mask when calculating the likelihoods. I n this anal- 
ysis we use the standard WMAP 'P06' mask (|Page et al.l 
|2003), which masks 26% of the sky. 



3 SIMULATED MAPS 

We generate simulated maps at the seven polarized nom- 
inal frequency channels for Planck (30, 44, 70, 100, 143, 
217, and 353 GHz). In our analysis, we do not apply 
beams or smoothing to the data; these would be included 
in a more realistic analysis but are not expected to signif- 
icantly affect results. Realizations of the CMB are gener- 
ated from a power spe ctrum computed usin g ACDM cos- 
mological parameters (|Komatsu et al. 1 120111 ). with either 
r = 0orr = 0.1. Diagonal white noise realizations are 
generated based on the noise levels ta ken from the Planck 
Bluebook ()Planck Collaboration! |2006| ). and we scale the 
given noise levels at beam-sized pixels to the correspond- 
ing noise level at A'aide = 16 sized pixels, with side 3.6°. 
This noise model is over-simplified as it contains no 1//- 
noise or other spatial correlations that are reported in the 
'early' Planck papers, which would increase effective no ise 
levels (jPlanck HFI Core Tearnll201ll : IZacchei et al. II2OIII '). 

For the foreground components, we use two baseline 
tests to benchmark the level of bias in the mismatch tests. 

In Test 1 (baseline with uniform /?s), spectral indices 
given by simple power-laws are used to simulate the syn- 
chrotron and dust foregrounds, and as a model in the compo- 
nent estimation. The simulated synchrotron Q and U emis- 
sion maps are modelled as power-law and given as an ex- 
trapolation in frequency of the polarized 23 GHz WMAP 
map: 



MP)=Q.siP){^y'''' 






(7) 



(8) 



We set the synchrotron spectral index to Ps = —3 uni- 
formly over the whole sky, consistent with o bservations 
by WMAP (|Page et all l2007l : ICold et al.|[2009h . The sim- 
ulated thermal dust Q and U emission maps are also mod- 
elled as power-law emission and generated by extrapolat - 
ing the predicted 94 GHz map in iFinkbeiner at al.l l|l993 ) : 
5i/(p) = 594 (p) (p;) '' ■ To generate the dust polarization 
angles we use a software package called the Planck Sky 
Model (PSM, version 1.6.6) developed by the Planck Work- 
ing Group 2. They closely match the sychrotron angles. The 
dust polarization fraction is set at 12%, which is scaled 
by a geometric depolarization factor due to the expected 
magnetic field configuration, resulting in an observed po- 
larization fraction of ~ 4%. We set the dust spectral index 



to Pd ~ 1.5 uniformly over the whole sky. This is consis- 
tent with polarization observations by WMAP at frequen- 
cies below 100 GHz, although at higher frequencies ther- 
mal emission is observed to devi ate from power-law (e.g., 
iPlanck Collaboration XXIVll201ll '). 

For the parametric model, we assume that the spectral 
index of the Galactic components do not vary over the fre- 
quency range considered, so the coefficients are given by 

0:2(2', /32) = dia,g[{iy/i/3o)'^^] 
0:3(2^, A) = dia,g[{iy/ 1/353)'^^]. 



(9) 
(10) 



Here we have defined the two spectral index vectors /32 and 
Pa for synchrotron and dust, respectively. We set the pivot 
frequencies to 30 GHz and 353 GHz. We impose Gaussian 
priors on the spectral index parameters of /32 = — 3.0±0.3 for 
synchrotron and P3 — 1.5 ± 0.5 for dust. The priors we have 
chosen have central value and standard deviation at approx- 
imately the average and range of values typic ally observed 
and pr edicted theoretically (s ee, for example, iFraisse et al.l 
(200^ : lDunklev et al.1 l|2009bl ') for further discussion). 

In Test 2 (baseline with non-uniform /?s), simple power- 
laws are again used to both simulate the foreground com- 
ponents and also as a model in the separation estimation, 
but the synchrotron index varies spatially over the sky. Dust 
emission is simulated as in baseline Test 1 but synchrotron 
emission is modelled as power-law with a spatially vary- 
ing Ps- The degree of spatial variation in the polarization 
spectral index has not yet been wel l-measured, but a realis- 
tic mo del is taken to be model 4 of iMiville-Deschenes et al.l 
(i2008l ). given by 



Ps 



log{P23/gfsS408) 



(11) 



log(23/0.408) 

where P23 is the WMAP polarization map at 23 GHz, g 
is a geometrical reduction factor (refiecting depolarization 
due to magnetic field structure), fs is the intrinsic polar- 
ization fraction from the cosmic ray energy spectrum, and 
5408 is the 408 MHz map of Haslam et al. (1982). The val- 
ues of Ps range from —3.3 to —2.8. The parametric model 
is as described in baseline Test 1, where we fit to power-law 
synchrotron and dust components. 



4 MISMATCH TESTS 

The set of tests described below are given a label identi- 
fier (A through I) and a short descriptive name to help 
the reader understand the results. In each test, we describe 
the model used to simulate the Galactic foreground compo- 
nent maps (known as the simulation) and then we describe 
the model used for the parametric component separation 
(known as the model). The mismatch tests are summarized 
in Table [1] We categorize the mismatch tests into the follow- 
ing three categories: incorrect model ( H4.ip : extra simulated 
components ( i]4.2|) : incorrect priors ( i]4.3|l . 

In every case, we define the parametric model for 
the sky signal using equation [3] Given that the CMB ra- 
diation is blackbody, the coefficient for cti is given by 
ai{i',Pi) — ai{i') = f{v)l, where the function f{iy) con- 
verts the CMB signal I from thermodynamic to antenna 
temperature. Though the spectral indices for Q and U in 
a given pixel are expected to be similar (following from the 
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Label 


Name 


Simulation 


Model 




Baseline Tests 


1 


Baseline uniform /3s 


sync power-law /3s = -3 
dust power-law /3d = 1.5 


sync power-law /3s = 
dust power-law /3d = 


-3 ±0.3 
1.5 ±0.5 


2 


Baseline non-uniform 0s 


sync power-law /3s = —3.3 to —2.8 
dust power-law 0^ = 1.5 


sync power-law /3s = 
dust power-law /3d = 


-3 ±0.3 
1.5 ±0.5 


Incorrect Model 


A 


Dust 2-component-a 


sync power-law /3s = -3 
2-component dust 


sync power-law /3s = 
dust power-law /3d = 


-3 ±0.3 

1.5 ±0.5 


B 


Dust 2-component-b 


sync power-law /3s = -3 
2-component dust 


sync power-law /3s = 
1-component dust 


-3 ±0.3 


C 


Synchrotron curvature 


sync curvature 

dust power-law 0^ = 1.5 


sync power-law /3s = 
dust power-law /3d = 


-3 ±0.3 

1.5 ±0.5 


Extra Components 


D 


1% Free- free 


sync power-law /3s = -3 
dust power-law /3d = 1.5 
1% polarized free-free 


sync power-law /3s = 
dust power-law /3d = 
no free-free 


-3 ±0.3 
1.5 ±0.5 


E 


1% Spinning dust 


sync power-law /3s = -3 
dust power-law /3d = 1.5 
1% polarized spinning dust 


sync power-law /3s = 
dust power-law /3d = 
no spinning dust 


-3 ±0.3 
1.5 ±0.5 


Incorrect Priors 


F 


Strong I3s prior mismatch 


sync power-law /3s = —3.3 to —2.8 
dust power-law /3d = 1.5 


sync power-law /Js = 
dust power-law /3d = 


-2.5 ±0.5 
1.5 ±0.5 


G 


Weak 13s prior mismatch 


sync power-law /3s = —3.3 to —2.8 
dust power-law /3d = 1.5 


sync power-law /3s = 
dust power-law /3d = 


-2.8 ±0.5 
1.5 ±0.5 


H 


Strong /3rf prior mismatch 


sync power-law /3s = —3 
dust power-law /3d = 1.5 


sync power-law /3s = 
dust power-law /3d = 


-3 ±0.3 
2.0 ±0.5 


I 


Weak /3d prior mismatch 


sync power-law /3s = —3 
dust power-law /3d = 1.5 


sync power-law /3s = 
dust power-law /3d = 


-3 ±0.3 
1.7 ±0.5 



Table 1. Summary of mismatch tests. 



assumption that the polarization angle does not change with 
frequency), unless otherwise stated, we allow the option for 
the indices to be sampled independently for Q and for U. 
Thus, our model is completely described by 6A'j, amplitude 
parameters A = (Ai,A2,A3) and iNp spectral index pa- 
rameters f3 = {l32 , fif , (32 , f^a ) ■ We impose a flat prior on 
amplitude-type parameters and Gaussian priors on the spec- 
tral index parameters. The model is estimated from UNp 
data points (seven frequencies with two Stokes parameters). 



We plot the likelihood curves for the estimated parame- 
ters, r and r, for each mismatch case and show the compari- 
son likelihood curve from its corresponding baseline test. By 
holding all parameters constant, except for the mismatch be- 
ing tested, we are able to quantify the level of bias induced 
by each type of mismatch. In this section we describe each 
test and present the numerical results; in Section [5] we dis- 
cuss their implications. 



4.1 Incorrect model 

Here we consider a subset of cases where the frequency de- 
pendence of the synchrotron and dust emission are modelled 
incorrectly. 



4-. 1.1 Thermal dust frequency dependence 

Thermal dust emission is well-approximated by a modi- 
fied black-body, with intensity scaling as u^B,y{T), where 
BuiT) is a black-body spectrum with temperature T. In the 
Rayleigh- Jeans limit, this approximates to the power-law as- 
sumed in our baseline simulations. Over a broader frequency 
range, the power-law approximation breaks down, and mod- 
elling the curvature becomes important. In the simplest ex- 
tension to a power law, it is common to fit for one or two 
parameters to describe the integrated dust emission from 
any line of sight: either the emissivity index /3, or emissivity 
plus temperature T. More realistically, the integrated dust 
emission arises from dust grains at various temperatures, 
so could b est be represented by th e sum of modified black- 
bodies. In iFinkbeiner at al.l (| 19991 ). a model with just two 
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Figure 1. Recovered distributions for the tensor-to-scalar ratio, r, for mismatched simulation and models, comparing the baseline (test 
1) with three mismatched cases for r = (left), and r = 0.1 (right). Modelling a two-component thermal dust simulation (modified 
black-body emission with dust at mean temperatures 16 K and 10 K) with a power-law dust spectral index (test A) biases r high by 
about Icr, as does neglecting a curvature in the synchrotron spectral index (test C). Modeling a two-component dust simulation with a 
one-component modified black-body model (Test B) has only a minor effect. 



components at mean temperatures 9.6 K and 16.4 K was 
found to be a good fit to the IRAS data. 

Here we consider two mismatches between model and 
simulation. In Test A (two-component-dust-a), the dust 
emission is simulated with two temperature components, 
while the para metric model fits t o a d ust power-law. We 
use model 7 of iFinkbeiner at all l| 19991 ). with [Q,U]{u) oc 
Aityl^^ B^{Ti) + A2u'^'^ B„{T2). In this model the first compo- 
nent is sub-dominant, with A2/A1 = 24.6. The dust emis- 
sivity indices are /3i — 1.5, /32 = 2.6 over the whole sky. 
Synchrotron emission is simulated as power-law with a spa- 
tially uniform jSa- The parametric model fits to power-law 
dust and synchrotron, neglecting the curvature of the dust 
spectrum. In Test B (two-component-dust-b), dust emission 
is again simulated with two temperature components (as in 
Test A), while the parametric model fits to a one-component 
dust model, [Q,U]{u) oc i'^B,j{T). We fix the temperature 
T over the sky to the values of T2 from the simulation, and 
estimate a single index Pd in every pixel. 

Using these test cases, we perform component separa- 
tion and use the resulting CMB maps to compute the like- 
lihoods for parameters t and r for the r = and r — OA 
simulations. The distributions are shown in Fig. [l] and re- 
covered mean values for r, and r, for these and all other 
tests are summarized in Table [21 For r = we quote 95% 
upper limits; for r = 0.1 and t we give 68% confidence lev- 
els. For r = 0.1 we find a non- negligible bias on r of la 
high for Test A, fitting a two-component dust model with a 
power-law, and a similar bias high for the optical depth, r. 
Using a one-temperature component model to fit the two- 
component simulation (Test B), recovers r with only ~ 0.2a 
bias. We see a similar effect for the r = case, where for 
Test A the recovered r value is greater than zero at Icr, but 
Test B is consistent with the baseline case. 



4- 1-2 Synchrotron frequency dependence 

Synchrotron emissi on is expected to be roughly power-law in 
frequency (see e.g.. lRvbicki fc Lightmanlll979l ). the result of 



relativistic cos mic-ray electrons a ccelerated in the Galactic 
magnetic field (jStrong et al.ll2007l ). However, a steepening of 
tlic index with frequency is also exp ected, due t o increased 
energy loss of the electrons (e.g., iBandav fc W olfcndale I 
I1991I : IStronget al.ll2007l ). The WMAP data are consistent 
with power-law emission, but a modest steepening would fit 
the data, and can be parameterized by a curvature of the 
spectral index. In a pessismistic scenario, the degree of steep- 
ening could vary significantly over the sky, or the frequency 
dependence could be ill-fit by a single curvature parameter. 
In Test C (synchrotron curvature), the simulated Galac- 
tic foreground includes a steepening of the synchrotron index 
with frequency while the parametric model retains power- 
law synchrotron emission. The synchrotron emission has 
spectral curvature such that the index decreases by 0.3 above 
23 GHz. Figure \T\ and Table [5] show the results from this 
third test case. The effect on the recovered CMB is non- 
negligible. We find that a synchrotron curvature simulation 
generates a bias of about la high in r, or 5r ~ 0.03, roughly 
the same level as the two-component dust simulation with 
power-law model. This mismatch also results in a 1.5a pref- 
erence for r > for the r = model. 



4.2 Additional polarized components 

Our model and simulations contain only synchrotron and 
thermal dust emission components. Other emission com- 
pon ents are not e x pecte d to be significantly polarized (see 
e.g.. iFraisse et al.l (|2008l ). and Section [5] for further discus- 
sion). However, both free- free and spinning dust emission 
are detected in intensity, and t hey may be mi n imally po- 
larized at the few-percent level. iMacellari et al.l l|201ll ) find 
an upper limit on spinning dust of 5% and an upp e r limi t 
on free-free polarization of < 3%. [Dickinson et al.l ([201l|); 
iLopez-Caraballo et al.l (|201ll ) reduce the upper limits on 
spinning dust polarization to ~ 1 — 2%. 

Test D (free-free) simulates a Galactic foreground that 
includes a 1% polarized free-free emission in addition to 
the synchrotron and dust emission. Free-free Q and U emis- 
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Test 


Recovered r Recovered r 
r = r = 0.1 


Bias (a) 


Recovered r 
r = 0.1 


Bias (ct) 


Baseline Tests 


(1) Baseline (uniform /3a) 

(2) Baseline (non-uniform 0^) 


< 0.03t 0.092 ± 0.033 

< 0.03t 0.092 ± 0.033 


— 


0.094 ± 0.005 
0.094 ±0.005 


- 


Incorrect Model 


(A) Dust 2-component-a 

(B) Dust 2-component-b 

(C) Synchrotron curvature 


0.02 ±0.016 0.125 ±0.037 

< 0.04t 0.096 ± 0.036 

0.03 ±0.020 0.125 ±0.039 


+0.9 
+0.2 
+0.9 


0.097 ± 0.005 
0.094 ± 0.005 
0.097 ±0.005 


+0.6 

<+0.1 

+0.6 


Extra Components 


(D) 1% free free 

(E) 1% spinning dust 


< 0.03t 0.091 ± 0.032 
<0.04t 0.094 ±0.033 


< -0.03 

< +0.03 


0.094 ± 0.005 
0.094 ± 0.005 


< +0.1 

< -0.1 


Incorrect Priors 



(F) Strong 0s prior mismatch 

(G) Weak 0s prior mismatch 
(H) Strong 0^ prior mismatch 
(I) Weak /3d prior mismatch 



0.168 + 0.047 0.197 ±0.047 

0.029 ±0.021 0.117 ±0.039 

0.133 ±0.044 0.224 ±0.040 

< 0.04t 0.111 ±0.034 



+2.1 


0.104 + 0.006 


+1.7 


+0.6 


0.096 ± 0.005 


+0.4 


+3.3 


0.107 ±0.005 


+2.6 


+0.6 


0.096 ± 0.005 


+0.4 



Table 2. Marginalized estimates and corresponding biases for r for simulations with 
T = 0.1. f These values are the upper 95% confidence levels for r = 0. 
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Figure 2. Recovered distributions for the tensor-to-scalar ratio, r, for simulations containing polarized components that are neglected 
in the models. The baseline results (test 1) are compared to those with a 1% polarized free-free component (test D), and a 1% polarized 
spinning dust component (test E), for r = (left), and r = 0.1 (right). At this polarization level, these components are sufficiently 
sub-dominant that they do not bias the recovered parameters. 



sion are given by Qff{u) — 0.01/ff (i^) cos(27) and C/ff(i^) — 
0.017fr(t^) sin(27), where Isiv) is a free-free intensity map at 
frequency v and 7 are the thermal dust angles. This assumes 
that the free-free polarization angles match the thermal dust 
angles, which is unrealistic but should not significantly af- 
fect conclusions. The free-free intensity is generated from 
the PSM, which is consistent with WMAP data. The para- 
metric model fits for power-law synchrotron and dust but 
omits the free- free component. 

Test E (spinning dust) includes a 1% polarized spinning 
dust emission in addition to synchrotron and thermal dust. 
Spinning dust Q and U emission are given by Qsdii^) = 
0.01/ad(j^) cos(27) and Usd{i^) = O.Ol/sdCi^) sin(27), where 
Isd{i^) is a spinning dust intensity map at frequency 1/ es- 



timated from the PSM, and the angles 7 are the same as 
the thermal dust angles. The parametric model omits the 
spinning dust component. 

The resulting likelihoods are shown in Fig. (2] and pa- 
rameters given in Table [21 We find that these small unmod- 
elled components have a negligible effect on the estimated 
parameters; the induced biases are within 0.04(7 of the base- 
line measurement in each case. 



4.3 Incorrect priors 

In our baseline model estimation we imposed Gaussian pri- 
ors of j3s = —3 ± 0.3 for the synchrotron spectral index, and 
Pd ~ 1.5 ± 0.5 for the thermal dust emissivity index. This 
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Figure 3. Recovered distributions for r, if prior distributions are imposed on spectral indices that do not exactly match the simulation 
inputs. The baseline (test 2) has a dust index /3^ = 1.5, and a synchrotron index with mean /3s = —3 over the sky. Indices for Stokes Q 
and U are fit in 3-degree pixels over the sky, with Gaussian priors /3s = —3 it 0.3 and [3^ = 1.5 ± 0.5. Offsetting the synchotron prior by 
1(7 to —2.5 ± 0.5 (test F), significantly biases the recovered r high (top panels, for r = 0, left, and r = 0.1, right). A ~ 0.5cr offset (test 
G) results in a smaller but non-negligible bias. Similar biases are found for offsets in the dust prior (bottom), for /J^ = 2.0 ± 0.5 (test H) 
and /Sfj = 1.7 ± 0.5 (test I). These biases arise from over-parameterizing the model in low signal-to-noise regions. 



allowed an estimate of the CMB in areas of the sky with a 
low signal-to-noise ratio. Even with seven frequencies, if the 
signal-to-noise ratio is low, the synchrotron and dust com- 
ponent can become degenerate with the CMB unless priors 
are imposed. 

The priors are astrophysically motivated; synchrotron 
emission is expected to have an index in the typical range 
—3.5 ;$ /?s ^ —2.5, depending on the injection spectrum and 
nature of diff usion and cooling (|Rvbicki fc Lightmanlll979l : 
iFraisse et al.ll2008i ). Thermal dust emission is expected to 
have emissivity in dex in the range 1 ;S /3 ^ 2.5 (see e.g., 
iFraisse et al.ll2008l ) . The 2a range of the prior therefore cap- 
tures physically reasonable beheaviour. However, our sim- 
ulations are perfectly matched to these priors: the simu- 
lated synchrotron indices are either exactly —3.0 in Test 
1, or have a mean over the sky of —3 in Test 2, and the 
dust was simulated to have an index of 1.5. The real sky 
will likely not match so well: we expect the emission to lie 
in the prior r ange, but will not precisely match the mean. 
iDickinson et a l. (2009) conducted a similar study to quan- 
tify the effect of priors using real data. Though they found 
that the priors had a small impact on the CMB spectra, they 
considered unpolarized emission, where foregrounds are rel- 
atively smaller. 



We test the effects of these prior choices by fixing the 
simulation spectral behavior, but choosing alternative Gaus- 
sian priors with means that are offset from the simulation 
inputs. 

Test F ('strong' /3s prior mismatch) examines a reason- 
ably strong case of mismatch between the model prior and 
simulation for synchrotron. Using Test 2 as the baseline, it 
simulates synchrotron emission with values of /3s that range 
between —3.3 and —2.8, but the parametric model assumes 
power-law synchrotron with a prior on /3s of —2.5 ±0.5. Test 
G ('weak' /3s prior mismatch) assumes a prior of —2.8 ±0.5. 
Test H (strong /3d prior mismatch) has a mismatch between 
the model prior and simulation for dust. Using the base- 
line simulations, the dust emission has /3d — 1.5 while the 
parametric model assumes a prior on /?<; of 2.0 ± 0.5. Test I 
('weak' /3d prior mismatch) assumes a prior of 1.7 ± 0.5. 

The likelihoods for these cases are plotted in Fig. |3l 
with parameters reported in Table [S] These mismatches 
result in the most significant biases. For synchrotron, the 
strong mismatch case results in a 3.5o" spurious detection of 
r (0.17 ± 0.05), for a model with no tensor component. The 
recovered value for r is also biased about 2cr high for the 
r = 0.1 case, and the optical depth r is high by almost 2a. 
The weak mismatch case, with prior —2.8 ±0.5, is biased by 
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~ 0.6(7 in r, with a spurious signal at the la level. Similar 
results are seen for the dust emission. For the strong mis- 
match a signal is significantly detected at 3(t when r = 0, 
and biased more than 3(t for r = 0.1 (returning 0.22±0.04). 
The weak mismatch case suffers from a bias of 0.6o- in r, 
and 0.4cr in r. 



5 DISCUSSION 

We have found that modelling polarized Galactic fore- 
grounds incorrectly can lead to significant biases in the re- 
covered CMB signal. In this section we discuss the reasons 
these biases are observed, and how they might be mitigated. 



5.1 Effect of priors 

When marginalizing over foreground uncertainty using a pa- 
rameterized method, components are distinguished by their 
frequency dependence. This provides a way of separating the 
black-body CMB signal from the foreground components. In 
the low signal-to-noise regime a prior on this spectral behav- 
ior breaks the degeneracy between CMB and foregrounds. 

However, we find that choosing an incorrect, yet phys- 
ically reasonable, prior for the frequency dependence can 
have a significant impact on the estimated cosmological sig- 
nal. With a simulated synchrotron spectral index between 
—3.3 and —2.8, and a Gaussian prior of —2.5 ±0.5 on the in- 
dex in each pixel, the tensor-to-scalar ratio is overestimated 
by ~ 3cr for an r = 0.1 model, or a spurious detection made 
when r = Q. The effect is less extreme when the mean of the 
Gaussian prior is closer to the input, —2.8, but a bias of la 
is still observed. In the limit of a low signal-to-noise ratio, 
this can be understood as equivalent to setting the spec- 
tral index to the wrong value over the whole sky. A prior 
of /3s = — 2.5±0.5 results in an index that is everywhere 
~ —2.5, instead of the mean simulated value 13 s ~ —3. Simi- 
larly, a prior on the dust index, or emissivity, oi Pd = 2.0±0.5 
results in an index of ~ 2.0 instead of the simulated 1.5. 

This incorrect recovery in regions having a low signal- 
to-noise ratio is demonstrated in the left panels of Fig. |4] 
for the synchrotron Q-Stokes component. Away from the 
Galactic plane, the index is estimated to be roughly —2.5 ± 
0.5. We also show in Fig. [5] the frequency dependence of 
the components, rms averaged over the masked sky in 3.6° 
pixels, and compared to the CMB signal in both E-modes 
and B-modes for r = 0.1. Assuming that the synchrotron 
pivot is fixed at 30 GHz, an index that is too shallow by 
/3s ~ 0.5 overestimates the synchrotron power by of order 
0.1 fiK in antenna temperature at the foreground minimum 
of 100 GHz. This is significant compared to the r = 0.1 B- 
mode signal, so a bias is expected. Similarly for dust, with 
a pivot at 353 GHz, a dust emissivity index too steep by 
Pd ~ 0.5 would underestimate the dust at 100 GHz by up 
to ^ 0.1 ^K in antenna temperature; significant compared 
to the r = 0.1 signal. 

This specific case where the prior is systematically dif- 
ferent to the input by up to Ict everywhere on the sky is 
a pessimistic scenario, but not implausible. To avoid the 
risk of bias, one must therefore take care in how the fore- 
ground model is parameterized. In the Bayesian framework, 
our chosen model has too many free parameters, given the 



low signal-to-noise ratio, so the result is being driven by the 
prior. To mitigate this, there are several ways of increas- 
ing the signal-to-noise ratio in the indices: including an- 
cillary data f rom complementa ry experiments like WMAP 
and C-BASS (IKing et al.ll2010l ). assuming common temper- 
ature and polarization spectral indices, using larger pixels to 
define the indices, or defining spectral indices in harmonic 
space to allow spatial coherence. 

We consider two of these possible improvements. Each 
three-degree pixel can have a distinct spectral index for I, Q, 
and U. The first natural improvement is to fix the Q and U 
spectral indices to be common in each pixel, /3q = PI/- Phys- 
ically this is reasonable; the polarized signal comes from the 
same region of the Galaxy for both Q and U-type, and can 
be expected to have the same frequency dependence, con- 
sistent with obs ervations (jKogut et al.ll2007l : iDunklev et al.l 
l2009al : lGold et al. 2009). We repeat Tests F and G with this 
condition (Tests F2 and G2), and show the recovered in- 
dex map in Fig. [l] with the likelihoods for r in Fig. |6] The 
index map now has a higher signal-to-noise ratio, and the 
bias on r reduced from more than 2a to la (for a prior of 
/3s = —2.5 ± 0.5). Fixing the temperature and polarization 
indices to be common is less physically motivated so we do 
not consider this here; depolarization effects could lead to 
different regions of the Galaxy contributing to the integrated 
polarization signal. 

The signal-to-noise ratio can also be improved by 
adding ancillary data that better traces the foregrounds. 
Since the synchrotron signal dominates at lower frequencies, 
additional data at the low frequency range will increase the 
synchrotron signal-to-noise ratio. We repeat Test F again 
(F3), adding simulated data from the WMAP 23 GHz K- 
Band channel, and projected C-BASS data at 5 GHz, to the 
simulated PZancfc data from 30-353 GHz. Figure |4] shows the 
significantly improved estimate of the synchrotron index in 
this case, which translates into a reduction in bias on r from 
2(T to la for a prior of /3s = —2.5 ± 0.5. With the low fre- 
quency data, the indices are better constrained by the data. 

A final obvious way to reduce the model freedom is to al- 
low less spatial variation in the indices. In the lim it of no spa- 
tial variation, this reduc e s to template cleaning (|Page et al.l 
|2007| ; iKogut et all l2007l : lEfstathiou et al.ll2009l ). with one 
spectral index over the whole sky. However, a concern with 
these methods is that they may not capture realistic spatial 
variation. The optimal balance is likely in between, requir- 
ing fewer than ~3000 parameters to describe the spatially 
varying frequency dependence. Such an app roach has been 
considered for polarization analysis in e.g., iDunklev et al.l 
pOOQah . where 48 synchrotron spectral index parameters 
were used for WMAP component separation. In making this 
choice with real data, it will be important to test that results 
do not depend on the prior placed on frequency dependence. 
If so, the number of parameters should be reduced, or ex- 
ternal data included where available. 



5.2 Effect of over-simplified model 



In Section r4. II we found that over-simplifying the frequency 
dependence of the two components can also lead to a bias in 
recovered parameters. Modeling the synchrotron as a power 
law everywhere on the sky, when it actually has a spectral 
curvature of C = —0.3, results in a ^ 0.03 bias high in r. 
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Figure 4. Estimated synchrotron spectral index for the Q-Stokes parameter (showing mean, top, and uncertainty, bottom), for a 
simulation with mean /3s = —3 and prior —2.5 ± 0.5 (Test F). Allowing free Q and U spectral indices, and using just 30-353 GHz 
data (left), the prior of —2.5 is returned in low signal-to-noise regions. If Q and U signals are assigned a common index (centre), the 
signal-to-noise is increased. If low- frequency simulated data from WMAP (23 GHz) and C-BASS (5 GHZ) is added (right), the spectral 
index map is recovered with high signal-to-noise. 



a. 



10.00 


1 1 1 1 1 1 


: 




/1^c■^,./Tn^-3 


- 




[ ' ) ^mV^/ ->^) 




S. (F) S,,{u/30r^l 






- X.. (G) Sju/50r'-' 


_ 




vw* . ^_ . ^_ . ^_ . ^_ . ^_ . . If 1 ^^ \ / K~^ j^^ ( 1 1 i^ \ y j^^ T 1 i v^ /m 






^^■. \^ J -^yiic L^urvaiure 




^'-■■- f ^ ) n i't- ''^^-^V-^ 


^^ 




^•■. l U Ujjj^y/JoJj 




■^^••. (u\ n ii/ /^^^V-° ^^^ 




1.00 


s\s.'-. \" J L/353l^l^/ OO J^ ^^<rf^ 


r ■^■■■■- (') D353(i^/353)'^^^<<r^- 


- 




'Os' \ ••. ^^^^^■'' 


- 




"^ r\j^ ' ■ ^.'■-'-'''''^^^ ■■' 


- 










^-^^ ^^^T^ . 












"^ ^\ ^^^^^''^^ ^ .-■' 






- V-;--^'-^ "^ - "^ ■.,.■■ 


- 




^,,„-'--''^ ^- '^\~"n^'^~~~^*>-^ 






^""""^^^ ^ ^ ^-.'^V ^ ^^^^^^^^^ 










0.10 


^.-'-''''''''^ -^ '' -•* \ ^S. ^ ■ ^^"^""---v. 


— 










^.-'''''''''^ -^ "^ ■' ^- ^S. "^ '■ ^Xs. 


- 










- ^--'••'''''''^ ^ '' -■■' N \v '^ "-. ^\ 


_ 
















-^ "^ ••' ^ ^v "^ ■ ^v 






"" .-■■'" ^s \C"- ■•••■ 


\ 




-■' "^ x,^^ \ '•-_ 


\ 




- ■' "\ x,^^ \ 




0.01 


1 1 1 1 1 1 1 N 1 \ s 1 





100 
Frequency (GHz) 



Figure 5. Frequency scaling of the foreground components in the baseline simulation (test 1), rms averaged over the unmasked sky in 
3.7° pixels ((. ~ 50 scales), and compared to the CMB E-mode signal for t = 0.1 (solid blue curve) and B-mode signal for r = 0.1 (dashed 
blue curve). If an incorrect spectral index in synchrotron or thermal dust is assumed (e.g., by imposing a prior: tests F, G, H, and I), or 
a synchrotron curvature neglected (test C), the over- or undcr-subtraction of foregrounds at ~ 100 — 150 GHz is significant compared to 
an r = 0.1 signal. 
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As in Sec 15.11 this can be understood as an overestimation 
of synchrotron at the 100 GHz range by up to ~ 0.05/iK 
in antenna temperature, illustrated in Fig. [S] Since some 
steepening is expected from synchrotron cooling, a strategy 
to prevent this bias would be to additionally marginalize 
over a curvature parameter. If the estimated CMB power 
does not change significantly with its inclusion, and the cur- 
vature is consistent with zero, this would justify neglecting 
the additional complexity. 

While we have examined only the case for synchrotron 
having a negative spectral curvature, there is some evidence 
to suggest that the spectral curvature could be positive 
(e.g., Dickinson et al. 2009; do Olivcira-Costa ct al. 200a; 
iKogut et al.l 120071 ). This is not unexpected since multiple 
spectral components can give a flattening of the effective 
synchrotron index. With real data, a positive curvature as 
large as 0.3 could be realistically considered. 

At the high frequency end, thermal dust emission is 
typically modelled as a modified black-body, characterized 
by an emissivity and temperature, with I{u) oc v^B^{T), 
and similarly for Q and U. This corresponds to our 'one- 
component' dust model. A more complicated model has a 



sum of two or more components with different temperatures. 
In Sec 14.11 we found that modelling a two-component dust 
model as a one-component dust model has only a small effect 
on the estimated CMB signal. This reflects that the sum of 
two modified black-bodies, one sub-dominant, scales with 
frequency similarly to a single black-body. 

A larger bias was found for a modified black-body mod- 
elled as a power law. In this case we find a la shift in recov- 
ered r, with the power-law model typically over-subtracting 
dust. The effect is similar to neglecting synchrotron curva- 
ture. While it is unlikely in practice that the dust would 
be modelled as a pure power-law, it is possible that one 
could make the wrong choice for the dust temperature. In 
these tests we fixed the temperature to the input value that 
was common over the whole sky, and varied just the emis- 
sivity in each pixel. To check for a possible bias with real 
data, one would ideally additionally fit for the dust tem- 
perature. Another approach to determine the dust temper- 
ature would be to use the temperature data, including the 
higher-frequency unpolarized channels of Planck (545 and 
857 GHz), and IRAS/DIRBE data up to ~ 3000 GHz. The 
dust temperature could then be assumed to be common for 
the polarization data. 



5.3 Effect of neglected components 

We find that neglecting sub-dominant polarized free-free 
and spinning dust components has a negligible effect on 
the results. This can be understood from Fig. [T] The sim- 
ulations include a 1% polarized signal, with the rms sig- 
nal of each component, averaged outside the Galactic mask, 
shown to be sub-dominant to an r = 0.1 signal in the 
range u > 100 GHz. The true polarization of these com- 
ponents is unknown, but is not expected to exceed this 
level. Observations of the Ophiuchi and Perseus cloud limit 
the polarization of spinning dust to be less than 2% at 
20-30 GHz (jPlanck CoUaboration XX II2OI1I ). and WMAP 
observations limit it to less than 1% over the whole sky. 
T hese levels are consistent with the spinning dust model 
byiDraine & Lazarian (1999). For this mask, spinning dust 
polarization has a slightly larger effect on r than free-free 
polarization. The spinning dust component is currently the 
most uncertain, so will be worth re- visiting with real data. 
There are fewer observational constraints on the po- 
larization of free-free emission. However, it should be in- 
trinsically unpolarized because the scattering directions are 
random. Secondary polarization can be generated at the 
edg es of bright free-free features from Tho mson scattering 
(Rv bicki fc Lightmanlll979l : iKeating et al.l lT998l. but lead- 
ing to less than 1% polarization at high Galactic latitudes. 
We have not therefore considered larger polarization levels. 
We have also not con sidered more exotic compon ents, such 
as a polarized 'Ha ze' IIDobler fc Fi nkbeiner 1 120071 ). or mag- 
netic dust models (jDraine fc Lazarian.. 1999i ). 



6 CONCLUSIONS 

Extracting robust estimates for the tensor-to-scalar ratio 
rely on modelling and subtracting polarized foregrounds. 
Since the polarized CMB signal is many times smaller than 
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the foreground emission, the need to get this right is par- 
ticularly acute. Many methods have been considered and 
implemented for foreground removal, but given the lack of 
data, the simulations are usually simple in form. 

In this paper we have begun to quantify the impact on 
estimates of r of incorrect foregound modelling. The tests 
were aimed at a detection of a signal with r — 0.1, but the 
goal of future missions is to reach r = 0.01 or lower, so we 
also consider an r = model. We conclude that neglect- 
ing a non-power-law frequency dependence of foregrounds 
may have a non-negligible effect on r; whereas neglecting a 
small free-free or spinning dust component is likely not to. 
We found that over-parameterizing the spectral indices had 
significant consequences; in the limit of a low signal-to-noise 
ratio the result can be highly prior-dependent. 

We discussed methods of mitigating possible bias, 
through model comparison as more complexity is added to 
the foreground model, and through increasing the signal-to- 
noise ratio on spectral parameters by reducing their num- 
ber and using ancillary data. We did not cover all scenarios 
of mismatch, but the approach of checking the goodness- 
of-fit through model comparison, and checking for a de- 
pendence of results on priors should be generally applica- 
ble. We did not explore the effects of different masks al- 
thou gh this will be import ant to investigate with data (see 
e.g., iDickinson at alj 12009 1. Data from Planck and ground- 
based and balloon experiments will further elucidate the 
nature of the polarized foregounds and allow their mod- 
elling to be refined. For full-sky data from f uture ultra- 
high sensitivi ty experiments such as CMBp ol (jBock et al.l 
|2009|), COrEi'The COrE Collaborationll201ll ). and LiteBird 
i Hazumi et al.ll2008l 'l. the effects studied here will be more 
important as we push towards r = 10~^ — 10~^ levels. 
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